function diff = model_1(theta,param,aux)


    global  search  
    theta=abs(theta)/(1+abs(theta));
   
    %% Calculate costants for the natural wastage region 
    region='N';

    %% Calculate roots in the natural wastage region
    a=-1/2.*param.sigma^2; b= -(param.mu-1/2.*param.sigma^2+(1-param.alpha).*(param.s.*param.lambda+param.zeta));
    c=+param.r+param.s.*param.lambda+param.zeta;
    param.gamma_N_1 = (-b+(b.^2-4.*a.*c).^(1/2))./2./a;
    param.gamma_N_2 = (-b-(b.^2-4.*a.*c).^(1/2))./2./a;
    clear a b c
     
    %% Create bounds for G=m_u/m_l
    G_l=1+param.c.*(param.r+param.s.*param.lambda+param.zeta)./(param.omega_0);
    G_u=(1-1./param.gamma_N_1)./(1-1./param.gamma_N_2).*(1+param.c.*(param.r+param.s.*param.lambda+param.zeta)./(param.omega_0));

    tmp_G= fzero(@ (G) (param.omega_0 + rho_fun( 0, region, param).*param.c)./param.omega_0.*phi_fun( G, region, param) ...
                -G.*phi_fun( G.^(-1), region, param), [G_l,G_u], optimset('display','off'));

    if param.K == 0   % OJS model
        
        %% Natural wastage region
        % Solve for the boundaries and constants in natural wastage region
        
        param.G=tmp_G;
        param.m_l = param.omega_0./(1-param.omega_1)./rho_fun(0,'N',param)./phi_fun( param.G, 'N',param);
        param.m_h = (param.omega_0+rho_fun(0,'N',param).*param.c)./(1-param.omega_1)./rho_fun(0,'N',param)./phi_fun( param.G.^(-1), 'N', param);

        param.J_1 = - (1-param.omega_1).*theta_fun( param.G, 'N', param).*param.m_l.^(1-param.gamma_N_1)./param.gamma_N_1./rho_fun(1,'N',param);
        param.J_2 = - (1-param.omega_1).*(1-theta_fun( param.G, 'N', param)).*param.m_l.^(1-param.gamma_N_2)./param.gamma_N_2./rho_fun(1,'N',param);

        %% Expansion region
        x0= fzero(@ (x) delta_OJS(abs(x)+param.m_h,param), 10^(-2), optimset('display','off'));
        param.m_u=param.m_h+abs(x0);
        f = @(x) ( (1-param.alpha).*2./param.sigma^2.*param.zeta./x./q_0_fun(x,param));
        param.psi=(1./q_0_fun(param.m_u,param)+integral(f,param.m_l,param.m_u));
        param.mean_x = mean_fun(aux.n,'x',param);
        param.mean_lnw = mean_fun(aux.n,'lnwage',param);

        diff=param;
        
    elseif param.K>0 % Baseline model
        
        %% Natural wastage region 
        % Solve for the boundaries and constants

        % Inital guess of m_l which is restrict to be within the relevant range
        tmp_m_l = param.omega_0./(1-param.omega_1)./rho_fun(0,'N',param)./phi_fun( tmp_G, 'N', param);
        param.m_l = theta.*tmp_m_l;

        % Set constants J_1 and J_2 such that J_p and J=0
        param.J_2= param.m_l.^(-param.gamma_N_2).*(param.gamma_N_2./param.gamma_N_1-1).^(-1).*...
            ((1-param.omega_1).*param.m_l./rho_fun(1,'N',param) - param.omega_0./rho_fun(0,'N',param)- (1-param.omega_1).*param.m_l./rho_fun(1,'N',param)./param.gamma_N_1);

        param.J_1= param.m_l.^(-param.gamma_N_1).*(param.gamma_N_1./param.gamma_N_2-1).^(-1).*...
            ((1-param.omega_1).*param.m_l./rho_fun(1,'N',param) - param.omega_0./rho_fun(0,'N',param)- (1-param.omega_1).*param.m_l./rho_fun(1,'N',param)./param.gamma_N_2);

        % solve for m_h such that J=c            
        x0= fzero(@ (x) J_n(param.m_l+abs(x),param)-param.c, 0.1, optimset('display','off'));
        param.m_h = param.m_l+abs(x0);
        clear x0
        if J_n_p(param.m_h,param)<0 
            error('J\prime(m_h) negative check ml guess')
        end


        %% Replacement region R
        region = 'R';

        % Calculate roots in the replacement region
        a=-1/2.*param.sigma^2; 
        b= -(param.mu-1/2.*param.sigma^2);%+(1-param.alpha).*param.zeta
        c=+param.r;%+param.zeta
        param.gamma_R_1 = (-b+(b.^2-4.*a.*c).^(1/2))./2./a;
        param.gamma_R_2 = (-b-(b.^2-4.*a.*c).^(1/2))./2./a;
        clear a b c
        
        % boundaries        
        J=J_n(param.m_h,param);
        J_p=J_n_p(param.m_h,param);

        % Assign option constants in the replacement region    
        param.J_rep_2=param.m_h^(-param.gamma_R_2)./(1-param.gamma_R_2/param.gamma_R_1).*(J -((1-param.omega_1).*param.m_h./rho_fun(1,'R',param)-param.omega_0./rho_fun(0,'R',param)) ...
            -(param.m_h*J_p-((1-param.omega_1).*param.m_h./rho_fun(1,'R',param)-2.*param.c.*(1-param.alpha)./param.sigma^2.*(param.zeta+param.s*param.lambda)))/param.gamma_R_1);

        param.J_rep_1=param.m_h^(-param.gamma_R_1)./(1-param.gamma_R_1/param.gamma_R_2).*(J -((1-param.omega_1).*param.m_h./rho_fun(1,'R',param)-param.omega_0./rho_fun(0,'R',param)) ...
            -(param.m_h*J_p-((1-param.omega_1).*param.m_h./rho_fun(1,'R',param)-2.*param.c.*(1-param.alpha)./param.sigma^2.*(param.zeta+param.s*param.lambda)))/param.gamma_R_2);

        if search~=0.5
            % calculate constant for searcher distribution within replacement region
            if param.zeta > 0 && param.s > 0
                try
                    x0= fzero(@ (x) root_psi(x,1,param),[param.zeta./(param.zeta+param.s.*param.lambda)+0.001,0.9] , optimset('display','off'));
                catch
                    x0= fzero(@ (x) root_psi(x,1,param),[param.zeta./(param.zeta+param.s.*param.lambda)+0.001,0.8] , optimset('display','off'));
                end
            elseif param.zeta == 0
                x0=0.5;
                param.B=(2.*(1-param.alpha)./param.sigma^2.*param.s*param.lambda).^(-1)-log(param.m_h);
            elseif param.s == 0
                x0=0.5;
                param.B=1;
            end           

            % find top of replacement region m_e such that J'(m_e)=0
            param = root_psi(x0,0,param);
            param.psi = q_fun(param.m_l,param)./q_fun(param.m_u,param);    
        end
        
        if param.K>0 && search == 1            
            % Calculate difference beween J(m_e) and C+c
            diff = J_r(param.m_e,param)-(param.K+param.c);
        elseif search == 0.5
            diff = param.J_rep_2;
        else
            param = root_psi(x0,0,param);
            param.psi=q_fun(param.m_l,param)./q_fun(param.m_u,param);
            param.mean_x = mean_fun(aux.n,'x',param);
            param.mean_lnw = mean_fun(aux.n,'lnwage',param);
            diff = param;
        end
    end
end

